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ABSTRACT 

We apply methods from Baycsian infcrcncing and graph theory to a dataset of 102 mid- infrared 
spectra, and archival data from the optical to the millimeter, to construct an evolutionary paradigm 
for z < 0.4 infrared-luminous galaxies (ULIRGs). We propose that the ULIRG lifecycle consists of 
three phases. The first phase lasts from the initial encounter until approximately coalescence. It is 
characterized by homogeneous mid-IR spectral shapes, and IR emission mainly from star formation, 
with a contribution from an AGN in some cases. At the end of this phase, a ULIRG enters one of 
two evolutionary paths depending on the dynamics of the merger, the available quantities of gas, and 
the masses of the black holes in the progenitors. On one branch, the contributions from the starburst 
and the AGN to the total IR luminosity decline and increase respectively. The IR spectral shapes 
are heterogeneous, likely due to feedback from AGN-driven winds. Some objects go through a brief 
QSO phase at the end. On the other branch, the decline of the starburst relative to the AGN is less 
pronounced, and few or no objects go through a QSO phase. We show that the 11.2/im PAH feature is 
a remarkably good diagnostic of evolutionary phase, and identify six ULIRGs that may be archetypes 
of key stages in this lifecycle. 

Subject headings: methods: data analysis - methods: statistical - galaxies: evolution - galaxies: active 
- infrared: galaxies - galaxies: interactions 



1. INTRODUCTION 

Ultraluminous Infrared Galaxies (ULIRGS, objects 
with rest-frame 1-1000/im luminosities of > IO^^Lq) 
play a fundamental role in the cosmological evolu- 
tion of galaxies and large-scale structures. First dis- 
covered in signific ant numbers by t he Infrared Astro- 
nomical Satellite (ISoifer et al.l I1984D. they are almost 
invar iably me rgers (iFarrah et all l2001l : iBushouse et all 
I2002t iVeiUcu x et all 12002 . 200 6i), powered by star for- 
mation and AG N activity, with the star formation usu- 
ally dominating (iGenzel et allliggS": 'Veil leux et al.lll999l: 
RigQPOulou et al.lll999t Imanishi et al. 20071: IVega et all 



2008ir ULIRGs are rare at low redshift, with less than 
fifty at z < 0.1, but become much more numerous at high 
redshifts, reaching a density on the sky of several hun- 
dred per square degree at z >. 1 (|Rowa n- Robins on et al 



1997 



1993 



Dole et al.ll2"00ll: [Hughes et al.lll 998: Bargc r et al 



Borvs et al.ll2003i : iMortier et al.ll200a) . The high 



redshift ULIRGs appear similar in some ways to those in 
the local Universe, i n that many of th e m are starburst 
dominated mergers (iFarra h et al.l[200l : ICh apman et al.l 
2003; Small et al. 2004; Ta kata et al. 2006; Borvs et al. 
20061: iVahante et al.i2007ciBe rta et al. 2007: B ridge et all 
2007t iLonsda le et all l2009t~ H uang ct al.. ■2009t) . though 
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there are also signs of differences, e.g. a higher 
fraction of systems with no signs of interaction 
([Melbourne et al.ll2008D. systematic ally different mid-IR 



spectral shapes (IFarrah et al 
cal environments ("Blai n et al 
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200 



and overdense lo- 
IFarrah et al] l2006l : 

iMag liocchctti ct al. 2 007[l. Reviews of t h eir pr oper- 
ties can be found in ISanders fc Mirabel! ([19961 ) and 
ILonsdale etlll ([2006f ). 

The cosmological significance of ULIRGs makes a solid 
understanding of them at low redshifts important, but 
there remain several uncertainties over the life-cycle of 
low redshift ULIRGs. We know that the merger activity 
triggers star formation and AGN activity, but we do not 
know how long the starburst lasts, whether or not there 
exist distinct 'AGN plus Starburst' or 'AGN dominated' 
phases, the fraction of ULIRGs that become Quasars, or 
if there are multiple evolutionary paths that a ULIRG 
can take. 

In this paper, we explore a new approach to study the 
evolution of the low redshift ULIRG population. Since 
many of the best diagnostics of star formation and AGN 
activity lie in the mid-infrared, we start with a mid- 
infrared spectroscopic dataset of z < 0.4 ULIRGs, taken 
with the Infrared Spectrograph (IRS. lHouck et alll2004D 
on board the Spitzer Space Telescope ([Werner et al.l 
120041 : ISoifer"era l. 2008). To this dataset we apply two 
novel analysis methods to determine trends in mid-IR 
spectral shape across the sample. First, we present a 
Bayesian based estimator of the degree of similarity be- 
tween a pair of ULIRG spectra. This estimator uses data 
from every resolution element, marginalizing over mea- 
surement error, luminosity, and foreground obscuration, 
to produce Bayes factors that describe the degree of re- 
semblance between every possible pair of spectra. Sec- 
ond, we use methods developed using graph theory to 
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study interconnected groups of entities to produce a 'net- 
work' diagram that visualizes these Bayes factors across 
the whole sample simultaneously. We combine these re- 
sults with archival data to propose an evolutionary de- 
scription of ULIRGs in the low-redshift Universe. 

This paper is structured as follows. In fj2]wc describe 
the sample selection and data, outline the method by 
which we calculate the Bayes factors, and describe the 
network construction. In ^J3]we present the Bayes factors 
and network diagram. In S2]wc assess the robustness of 
the Bayes factors and the network diagram, and possible 
sources of error. Discussion can be found in S}5l We 
summarize our conclusions in fj6l We assume a spatially 
flat cosmology with Hq = 70 km s~^ Mpc~^, = 1, and 
= 0.3. 

2. METHOD 
2.1. The Sample 

The sample was observed as part of the IRS Guaran- 
teed Time program to obtain mid-infrared spectra of low 
redshift IR-luminous sources (Spitzer program ID 105), 
selec ted from the IRAS IJy (|Kim fc Sanders! |1998[ ) and 
2Jy IStrauss et al'lggO") surveys, and from the FIRST 
sample (Stanford et al. 2000i). The sample is slightly bi- 
ased towa rds sources with w arm infrared colors, as de- 
scribed in iDesai etahl (|2007[ ) , but should still be repre- 
sentative of the low-redshift IR-luminous galaxy popu- 
lation. A few of our sample have IR luminosities that 
fall outside the usual definition of a ULIRG, but for sim- 
plicity wc refer to all of our sample as ULIRGs for the 
remainder of this paper. 

Low resolution spectra (5.2^m - 38.5/im, R^^ 60 — 125) 
were obtained of 118 objects, and high resolution spec- 
tra (9.6/im - 38.0/im, R^ 600) were obtained of a subset 
of 53 objects. Data reduction methods a nd initial re- 
sults a re p resented in lArmus et al.l ()2004D ; ISpoon et al.l 
(|2004[ ) a ndlArmus et all (120061). Further results are pre- 
sent ed inlHigdon et all (|2006f l: ISpoon et all (|2006l . l2007t l 
and iDesai et al.l ()2007t ). Atlases of the high and low res- 
olution spectra can be found in lFarrah et al.l (|2007t ) and 
Armus et al (2009, in preparation), respectively. 

The high resolution spectra contain more elements 
than the low resolution spectra, but sky background can- 
not be subtracted from them as we lack dedicated con- 
temporaneous sky observations. Therefore, we use the 
low resolution spectra. We exclude objects with z > 0.4, 
so that we see approximately the same wavelength range 
for all objects. We also exclude IRAS 11119-^3257 and 
IRAS 233654-3604 as they have poor quality data in the 
Short-Low modules. This leaves 102 objects, listed in 
Table [TJ We smooth the spectra to the instrumental res- 
olution using a 0.4/zm boxcar, and assume a 5% flux error 
for each resulting resolution element, consistent with the 
observed variations between individual nod positions. 

2.2. Analysis 
2.2.1. Bayesian Measures of Resemblance 

To compute the level of resemblance between any pair 
of (rest-frame) spectra we adopt a Bayesian approach. 
For any two spectra, A and B, we compute: 

^ P(A,B|different) 

P(A,B|same) ^ ' 



where P(A, Bjsamc) is the probability density that the 
two spectra are identical, and P(A, Bjdifferent) is the 
probability density that the spectra are different . The 
quantity TZ is thus the Bayes factor'' (j Jeff revs! 119611 : 
IConnoUv et al.l l2006[ ) quantifying the belief^ that this 
pair of spectra arise from sources whose physical prop- 
erties (or at least those that give rise to the mid-IR 
emission) are the same. In essence, we are performing 
a pixel-by-pixel comparison between two spectra, where 
no spectral region is preferentially weighted. While this 
method could be used to compare data to models, we are 
here making no model comparisons, instead comparing 
the spectra to each other. 

The simplest use of this method would involve comput- 
ing TZ for every possible pair of 'raw' (i.e. reduced, rest- 
frame but otherwise unaltered) spectra. This, however, is 
not enough. There exist two variables that will increase 
the Ti. values, but which do not necessarily reflect physi- 
cal differences. The first is instrument noise; differences 
between spectra that arise due to Gaussian fluctuations 
in the measurements should not contribute to the TZ val- 
ues. The second is cold foreground extinction; if object 
A and object B are intrinsically identical, but object A 
has a thicker screen of cold dust in front of it, then this 
should not contribute to the TZ values either. We make a 
third, simplifying assumption that differences in mid-IR 
luminosity (i.e. multiplicative scalings between spectra) 
do not reflect 'real' differences. Therefore, in calculat- 
ing the TZ values, we marginalize with respect to instru- 
ment noise, mid-IR luminosity, and cold foreground ex- 
tinction. The resulting (fully marginalized) probabilities, 
P(A, Bjsame) and P(A, B|differen t), are thus measures 
of the evidence for the hypothesis (jSivial [19961 ). The full 
methodology is described in the Appendix. 

Finally, we adopt a boundary condition for log]^Q(7?,); 
pairs of spectra with log]^Q(7?.) below this boundary are 
treated as similar, and those pairs with log]^Q(7?.) above 
it are treated as different. We set this boundary at 
logig(7^) = 0. In frequentist ter ms, this boundar y is 
equivalent to demanding < 0.8 (jSellke e t af 20011). In 
Bayesian terms, using the scale given in [Jeffre ys (1961.) . 
this corresponds to 'marginal' strength of evidence. We 
explore the sensitivity of our results to the error in TZ 
and choice of boundary condition in fj4| 

2.2.2. Network Construction 

With a sample of 102 galaxies, we have ^"^(72 — 5151 
TZ values, one for each possible pair. Our second require- 
ment is a method to study these Bayes factors across the 
sample. 

This is an example of studying a pairwise-connected 
group of entiti es. Other example s include a computer 
network (e.g. ISiganos et all [2003[ ). predator-prey rela- 
tionships among animals, or 'social' networks such as 
friendships between individuals in a group. As such, a 
common terminology has arisen to describe them. Each 

^ Arguably, the ratio of posteriors - i.e. 
P(difTeront|A,B)/P(same|A,B) is a more intuitive statistic. 
However, the calculation of the odds requires one to make prior 
assumptions about P(same) and P(diffcrent), which we prefer to 
avoid. 

* We use the word 'belief in its Bayesian sense, i.e. the odds of 
a successful trial of the truth of a given proposition, and not in the 
colloquial sense 
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Fig. 1. — Six examples of the 'best' fits (for which —150 < log{TZ) < —250), calculated by minimizing Eqn. \Ml\ The solid (dashed) 
black line is the IRS spectrum for object A (B). The solid (dashed) red line is the predicted intrinsic spectrum (i.e. without foreground 
extinction) for object A (B). The solid (dashed) blue line is the predicted intrinsic spectrum with cold foreground extinction applied for 
object A (B). The identical shapes of the solid and dashed red line shows that the intrinsic spectra of objects A and B are the same (they 
differ by a factor of fe-f'^" /(I — /)e'^~^"='''"= ). The identical shapes of the solid (dashed) black and blue lines shows that the observed and 
'predicted observed' spectra of object A (B) are the same. We do not constrain the shape of the intrinsic spectrum, and are showing the 
purely mathematically 'best' fits, so some predicted features in the intrinsic spectra, such as the very strong silicate emission feature in the 
bottom left panel, may not be 'real'. 
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Fig. 2. — For comparison with Figure \u two examples of the 
'worst' fits calculated by minimizing Eqn7 lA3T1 where log{TZ) ~ 
3000. Here the the solid and dashed red lines are identical, but 
the blue and black lines are completely different for both objects, 
showing that the IRS spectrum and the 'predicted' IRS spectrum 
do not match each other. 



entity (e.g. a computer or a person) is a 'node'. Connec- 
tions between nodes arc 'edges' if they have no direction 
(for example, if two computers share data) and 'arcs' if 
they have direction (for example, lions eat gazelles but 
gazelles do not eat lions). The number of edges con- 
necting to a node is the 'degree' of that node. The 
nodes in our network are the ULIRGs. The pairs where 
log]^o(''^) < are connected via edges, while those pairs 
with logjQ(7?.) > are not connected. 

Several methods have been developed to plot nodes 
and their connections in informative ways. In our 
case we require a method that (a) places connected 
nodes close together, and (b) produces as few 'cross- 
ing' edges as possible. A suitable algorithm for this is 
a 'force-directed' algorithm, in which attractive and re- 
pulsive forces_govBra_tlwajra^ 

edges (|Kamada &: Kawailll989HFruchterman fc Reingoldl 
Il991t ). Connected nodes attract each other along edges, 
while all nodes repel each other. The attractive force 
is modeled as if the edges are springs (i.e. a Hookes 
law type force) while the repulsive force is modeled as if 
the nodes arc electrically charged (i.e. a Coulomb type 
force). The parameters of the two forces are adjusted, 
and nodes allowed to move according to the forces act- 



ing on them until (a) an equilibrium state is reached in 
which the positions of the nodes and edges do not change 
appreciably, and (b) the nodes and edges can be seen si- 
multaneously. 

To create the network for our sample we use two soft- 
ware packages; the Network Workbench tool^, and Cy- 
toscape^'^. 

3. RESULTS 

The pairs of sources for which log]^Q(7?,) < are given in 
Table[TJ Examples of pairs where logiQ(7?.) < are shown 
in Figure [ll and examples of pairs with logiQ(7?.) > are 
shown in Figure [21 Presenting all of the TZ values would 
take an unreasonable amount of space, so instead wc plot 
their histogram in Figure ID The adjacency matrix, A, 
for our network^ ^ is too large to give in full, but the first 
few elements of it read: 



A = 







1 

10 



where the diagonal elements arc zero as our network con- 
tains no loops. The network for the sample is shown in 
Figure [21 

3.1. Structure 

Figure [21 exhibits a strong degree of connectivity. Only 
four objects (16,41,48,79, which arc not plotted) are not 
connected to any other. All the other nodes are con- 
nected to at least one other node, with most having three 
or more connections. The average node degree is high, 
at 10.2, and the average shortest path between any two 
nodes is short for a 102 node network, at 3.2 edges. 

There is significant variation in node degree across the 
diagram (Figure [5]). We quantify this by computing the 
fe-Ne arest-Ncighbor distribution (jPastor-Satorras et al.l 
l2001h for our network. We find that, as the average de- 
gree per node, fcmean, increases, so does fcjvAr; k^j^j ~ 0.5 
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with the caveat of the small number of nodes. Figure [5l 
shows a correlation between the degree of a node and 
that of its neighbors - nodes with a high degree arc more 
likely to be connected to other nodes with a high degree 
- and is thus an 'assortative' network. 

We identify at least two substructures. The first is 
a strongly interconnected group centered on object 29 
(IRAS 03000-2919), accompanied by some outliers on 
the left hand side, containing ~60% of the sample. We 
call this group A. The second is a weakly interconnected 
group extending in the rightward direction from group 
A, and containing the remaining 40% of the sample. 
It is plausible, given the two 'branches' in this second 
group, that it is composed of two groups; one extend- 
ing along the top of Figure [3l and centered on object 15 
(IRAS 00275-2859), and one on the lower side of Figure 



sity 



This 
and 



tool is developed jointly by 
Northwestern U niversity, and i 



Indiana Univer- 
1 available from 



'http:/ /nwb.slis. In diana. edul 

Available from |http:/ /cytoscape.org/| 

The adjacency matrix, A for an undirected graph with n nodes 
is defined as the n X n matrix where Aij is the number of edges 
from vertex i to vertex j , and An is the number of loops for vertex 



Fig. 3. — The network for our sample (generated using a spring-embedded algorithm within Cytoscape). The numbered points, or 
'nodes', are the objects in Table [T] An edge between two nodes indicates that log{Ti) < for that pair of objects. 
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Fig. 4. — Histogram of the Bayes factors for every possible pair 
of spectra (5151 in total) in Table [T] computed using the miser 
algorithm. 

[31 centered on object 97 (Mrk 1014). We label these two 
branches groups B and C respectively. Solely with this 
diagram to go on however, this subdivision is tentative, 
as the purported Group C resembles the 'outhers' on the 
left side of group A. We discuss the robustness of this 
subdivision in ijS.ll 

3.2. Node properties 

In studies of networks, insight can be gained by cod- 
ing the nodes according to some p roperty of the nodes 
for ex ample, coding the networks in lLusseau fc NewmanI 
2004D by gender or age revealed clear sub-communities) . 
We adopt this approach in this section. We here present 
the individual coded diagrams, and interpret them in fJS] 
Optical Spectral Type: (Figure (5]) The majority 
of the objects in group A have HII or LINER optical 
spectra, with a few Sy2's. This pattern is reversed in 




Fig. 5. — Network diagram, with the nodes color-coded by the 
number of edges connecting them (black=<4, red=4-7, yellow=7- 
14, green=14-20, blue=>20). 



group B; most of the objects have Sy2 or Syl spectra 
with a small number of LINERS and HII's, especially 
towards the right hand side where nearly all the objects 
are Syl's. In group C there appear to be approximately 
equal numbers of all spectral types. 

IR Luminosity: (Figure [6]) We expect a weakened 
correlation with 1-1000/xm luminosity, given the large un- 
certainties caused by the paucity of flux measurements at 
> 100/im. This expectation appears to be borne out^^; 
low and high luminosity systems are found in all three 
groups in approximately equal numbers, and there are 
no clear trends. There are no high luminosity systems in 
group C, though this may be due to the small number of 

A weakening in the correlation should not arise from the 
marginilization over IR luminosity. Marginilizing over luminos- 
ity will remove any dependence on luminosity. However, if lumi- 
nosity itself depends on (say) spectral shape at long wavelengths 
then the dependence on luminosity will remain, since we have not 
margnilized over spectral shape at long wavelengths. 
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Fig. 6. — Network diagram, with the nodes color-coded by Left Panel: IR luminosity, and Right Panel: optical spectral class (black: 
unknown, blue: HII, green: LINER, yellow: Sy2, red: Syl) 
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Fig. 7. — Network diagram, with the nodes color-coded by Left Panel: projected nuclear separation (black: unknown, blue: >12kpc, 
green: 6-12Kpc, yellow: 0.1-6Kpc, red: single nucleus) Right Panel: black hole mass (black: unknown, blue: > 2.5 X 10*Mq, green: 
8.0 <Mx10^Mq < 25.0, yellow: 5.0 <MxlO''M0 < 8.0, red: < 5 X IO'^Mq) 



objects in this group. 

Projected Nuclear separation^'^: (Figure [7|) There 
are strong caveats in interpreting this network; the imag- 
ing is heterogeneous (e.g. ground-based for some, space- 
based for others), the separations are projected rather 
than real, premergers that are widely separated can 
be erroneously identified as single nucleus systems and 
vice versa, an d nuclear separations are degenerate with 
merg er stage (jBarnes fc Hernquistlll992l : iDubinski et al.l 
Il999t ). The figure does however show trends. While the 
single nucleus systems are found in all groups, the ma- 
jority of them are in group B, including nearly all the 
objects at the end. Conversely, the widely and moder- 
ately separated systems are found almost exclusively in 
group A; with a few in groups B and C. 

Black Hole Mass^^: (Figure [7|) Here we use only 
those BH masses measured via velocity dispersions (see 
e.g. iTremaine et al.l |2002[ ). but the caveats are even 
stronger than for the nuclear separations diagram. Only 
33 measurements are available, the random uncertainties 
are large, and the measurements depend critically on the 
calibration of the Mbh — M^uige relation. Most impor- 
tantly, the quantity that we're really interested in is not 
Mbhi but rather the increase in black hole mass dur- 
ing the merger (i.e. l^MBH\merger)- The instantaneous 



Taken from IRigopoulou et al.l [T999I: IFarrah et al I 
Meusinger et al. 2001: Cui et al. 2001; Bushous e et all 
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Veilleux et al . 2002, 2006,; Bianchi et al...2008 and rcscaled to our 



cosmology where appropriate 

1"' Taken from lDasvra et al.|[2006al lbl and IKawakatu et al.ll2007l 

and rcscaled to our cosmology where appropriate 



snapshots of black hole mass are of limited use since the 
black holes in the progenitors can in principle have a 
wide range of starting values. Therefore, this diagram 
is of limited interest. The intermediate mass black holes 
are spread randomly through the diagram, but 4/5 high 
mass black holes are at the end of group B, and 11/16 
low mass black holes are in group A, with the rest lying 
in the first part of group B, or group C. 

PAH Equivalent Width: (Figure H) The mid-IR 
spectra of many ULIRGs show broad emission features 
at 6.2/xm, 7.7/xm, 8.6//m, 11.2/im and 12.7/j,m, attributed 
to bending and stretching modes in neutral and ion- 
ized Polycyclic Aromatic Hydrocarbon (PAH) molecules, 
and it is now accepted that these features indicate ongo- 
ing star formation. Therefore, the prominence of PAH 
features above the continuum, which we quantify via 
equivalent width, is a crude but reliable measure of 
the energetic importance of star formation^^ (see also 
Genzel et al"l998; 'RigO DOulou et al.lll999t lArmus et all 
2006: Desai et al. 2007). As there is still debate over 
the use of individual PAH features as star formation rate 
diagnostics, we show networks coded by the 6.2/im and 
11.2/^m PAH features. The 11.2yLtm diagram is particu- 
larly striking. All the objects in group A have promi- 
nent PAHs and there is a high degree of homogeneity in 
their strengths. The PAHs then decline in prominence 
as we move left to right through groups B and C, until 
we reach the right hand side of both groups where the 

As opposed to PAH fluxes, which measure the absolute lumi- 
nosity of the starburst 
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PAHs are negligible. The 6.2/im diagram shows the same 
trends, but less obviously; most of the objects in group A 
have prominent PAHs with some outlying objects show- 
ing weakened features, and there is a less pronounced, 
though still clear decline in PAH strength as we move 
left to right through groups B and C. 

Sil icate strength: (F i gure 51 defined in ISpoon et al.l 
120071 and ISirockv et all boOsT The 9.7/im feature is 
thought to arise from large silicate dust grains; in absorp- 
tion when 'cold' silicate grains absorb mid-IR continuum 
emission from a background source, and in emission when 
the silicate grains are 'hot'. It is usually interpreted as a 
measure of the obscuration towards the central, sub-Kpc 
nuclear regions. Under this interpretation, a prominent 
silicate feature is more correlated with AGN activity than 
star fo rmation - an absorp tion feature suggests a buried 
AGN (jlmanishi et al.ll2007D and an emission feature sug- 
gests an unobscured AGN (jHao et al.l |2005[ ). There is 
a caveat in interpreting this network though - the sili- 
cate feature is broad, substantially more so than a PAH 
feature, so its contribution to the TZ values will be com- 
mensurately larger^^. We cannot reliably gauge the mag- 
nitude of this effect, so the conclusions that can be drawn 
from this diagram are limited. We do however see trends. 
The silicate strengths of group A are fairly homogeneous; 
nearly all are moderately to heavily obscured, with (per- 
haps) slightly higher values for the outliers. Group C 
and the first part of group B are more varied, with a 
wide range of silicate strengths. The nodes at the end 
of group B universally show negligible absorption, or sil- 
icates in emission. 

As with the IR luminosities, the marginilization over fore- 
ground extinction should not affect any correlations with silicate 
strength 




Silicate 
Strength 



Fig. 9. — Network diagram, with the nodes color-coded by silicate 
strength. 



4. RELIABILITY 

We assess the accuracy and precision, along with pos- 
sible sources of error, on both the Bayes factors and the 
network in this section. 

4.1. The Bayes Factors 

We use three methods to examine the behavior of the 
Bayes factors. First, we assess the precision and accu- 
racy of the whole procedure. We compare what happens 
to the logiQ{TZ) values in two situations; (1) when the in- 
put spectra are intrinsically different, and (2) when the 
input spectra are intrinsically identical but where one 
has been scaled in luminosity and/or dust extinction. 
In both situations we vary the bin-to-bin uncertainties 
from 2% to 50%. The results from these simulations 
are shown in Figure [TOl This shows that, as the un- 
certainties become smaller, there is a higher probability 
that (a) log]^Q(7?.) < for those ULIRGs whose intrin- 
sic spectra are the same, and (b) log]^Q(7?.) > for those 
ULIRGs whose intrinsic spectra are different. Hence, the 
Bayes factors are behaving as expected. When the un- 
certainties are large there is an almost 100% overlap of 
those ULIRGs with 'same' and 'different' intrinsic spec- 
tra. However, the distributions are not centered around 
zero. We have found that when the uncertainties are in- 
creased even further (~ 10000%), pairings of same and 
different intrinsic spectra have a mean logjQ(7?,) closer to 
1, indicating this effect is primarily due to an increasing 
lack of accuracy in the integration as the uncertainties 
become smaller and the parameter space becomes more 
sparse. 
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Fig. 10. — Histograms of the Baycs factors assuming five different Gaussian bin-to-bin uncertainties. The dark grey histograms are 
those spectral pairs for which the intrinsic spectra for the two ULIRGs are the same but the dust extinctions are different. The light 
grey histograms are those ULIRGs for which both the intrinsic spectra and dust extinctions are different. As the uncertainties become 
smaller, there is a higher probability that logjg(7^) < for those ULIRGs whose intrinsic spectra are the same, and that logjQ(7^) > 
for those ULIRG whose intrinsic spectra are different. Hence, the Bayes factors are behaving as expected. The increased separation 
between histograms for intrinsically different and intrinsically similar pairs as the errors decrease is somewhat obscured by the limited 
accuracy /precision of Ti from the Monte Carlo integration. 
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Second, we assess the accuracy and precision of the 
integration algorithm used to compute the Bayes fac- 
tors. We used the miser algorithm , which is (to ou r 
knowledge) the most robust available (jPress "erallfl99l . 
but as it is a Monte Carlo algorithm, it has uncertain- 
ties associated with it which are difficult to compute^^ . 
Therefore, we adopt a conservative approach. We com- 
pare in Figure [TT] the log]^Q(7?.) values obtained with the 
miser algorithm, and an alternative algorithm, VEGAS. 
The spread is large, with a la dispersion on logiQ{TZ) of 
~ 50. The actual uncertainties due to the limited accu- 
racy/precision of the Monte Carlo integration are, how- 
ever, much lower than this, as the VEGAS algorithm is 

Usually, one can approximate how close one is to the "true" 
integral through the variance of the integrand. However, in our 
case, the errors are ~ 100% because the parameter space is both 
sparse and dynamic. Further, it is known a priori that the inte- 
grand can never be less than 0. 



less robust than miser. Furthermore, the effect on Fig- 
ure [3] of uncertainties on individual TZ values is likely to 
be small, as we simply demand that log]^Q(7?.) < 0. 

Finally, we examine whether small regions of the spec- 
tra dominate the derived TZ values. This is important 
for assessing the reliability of some of the coded variants 
of Figure [Sj the PAH and silicate strength plots repeat 
(some) information already in the TZ values, and so trends 
may be artifically amplified. We selected 200 pairs at 
random, and removed from each spectrum a 0.7/im wide 
region centered on the 11.2/im PAH feature. While there 
was some variation in the derived TZ values, we found 
that, in 97% of the cases, \ogiQ{TZ) did not change sign. 
Wc repeated this test for 0.7/im wide regions centered on 
the 7.7/im PAH feature and in the continuum at 15/im, 
and found similar results. We conclude that, while spec- 
tral windows of width ^ 1/im do contribute to the TZ 
values, they do not dominate them. Therefore, we argue 
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that if a spectral feature spans a small wavelength range, 
and contains information on a specific physical property, 
then coding Figure [3] by that feature tells us how that 
property varies with network position, while introducing 
minimal contamination. 

4.2. The Network 

In this section, we assess four possible sources of error 
in the methods used to create Figure [3] 

First, Figure [3] could simply be randomly connected 
points, and not contain any 'real' structure. We assess 
this in two ways. Qualitatively, Figure [3] does not re- 
semble a random network with 102 nodes, examples of 
which are shown in Figure [T2l Quantitatively, it has been 
shown that randomly co nnected graphs have a Poisso- 
nian degree distribution (jErdos fc Renviiri959f ). The de- 
gree distribution for Figure [3] is shown in Figure [131 It 
is not well matched to a Poisson distribution. 

Second is the effect on Figure [3] from the limited preci- 
sion of the Bayes factors. The variation in the \ogi^{TV} 
values due to the Monte-Carlo integration is difficult to 
determine, though we showed in ij4.1l that the upper limit 
is ~ 50. So, to assess this, we show in the top panel of 
Figure [14] the layout obtained after randomly changing 
all the \ogiQ{TZ) values by ±50. Even when using the 
upper limit on the errors, we get essentially the same 
structure. We conclude that Figure [3] is robust to within 
the accuracy and precision of the TZ values. 

Third, one could argue that the structure of Figure [3] 
is an artifact of the priors, and does not refiect trends 
in the data. To examine this, we test the sensitivity of 
Figure [3] to the maximum allowed cold foreground dust 
extinction, which is the only prior we adopt. Figure [3] 
assumes a weak upper limit o n the cold foreground ex- 
tinction of Ag.j^ni ~ 50 (see ^A.4p . In the lower panel 
of Figure [14] we show the diagram obtained if this limit 
is relaxed further, to Ag.7^m — 80. It closely resembles 
Figure [3] though there are differences in the number of 
edges connecting some nodes. As we are using a Bayesian 
approach, the relaxation of the constraint on Ag y^m can 
cause the TZ values to rise or fall; for example, IRAS 
08572-1-3915 is now not connected to any other node and 
is not plotted. Overall, we conclude that changing the 
prior on Ag.y^m does not substantially affect the struc- 
ture of Figure [3] 

Fourth, it is possible that the algorithm used to gen- 
erate the network is a source of error. The methods 
described in §2.2.21 start each network from a random 
'seed' position, and assume parameters for the attrac- 
tive and repulsive forces. Therefore, the final appear- 
ance of a network will differ from case to case, and it 
may be that the seed position or the force parameters 
are governing the final appearance, not the input data. 
To test this, we checked to see if different seed positions 
or different force parameters gave networks without the 
structure seen in Figure [3] and found they did not. To 
illustrate, we show in Figure [T5l an example of a 'raw' net- 
work for our data, made using the default parameters for 
the spring-embedded algorithm within Cytoscape. The 
same structures can be seen in Figures 131 and [T51 the only 
difference being that the nodes in Figure [TSl are too small 
to identify by number. We conclude that our network is 
at least reasonably robust to the choice of parameters 
of the algorithm. This test illustrates a common prob- 
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Fig. 11. — A comparison of the log(7?.) values computed using 
two Monte-Carlo algorithms; miser (which we use in our analysis), 
and VEGAS. 

lem with this type of analysis; the 'raw' networks are 
not always amenable to visual interpretation. This is a 
particular problem in our case; we want all the nodes to 
be individually identifiable, but the node labels in the 
raw networks are invariably too small to see. Therefore, 
to arrive at Figure [3] we adjusted node/edge properties 
such that the information in the network was preserved, 
but in which the individual nodes can be identified. 

It is difficult to be certain of the robustness of Figure 
[3] as this technique has (to our knowledge) not been pre- 
viously used to examine any astronomical dataset, and 
so there exists no previous study for comparison. We 
have however comprehensively tested the robustness of 
Figure [3] and not found any significant problems. We 
conclude that this possibility is remote, and proceed on 
the assumption that the diagram reflects real trends in 
the data. 

5. DISCUSSION 

Figure [3] highlights intrinsic similarities in the mid-IR 
spectra of low redshift ULIRGs. It places little emphasis 
on individual spectral features (see ij4.ip . and we have 
marginalized over IR luminosity, foreground cold dust 
extinction and detector noise. Furthermore, as our sam- 
ple selection is almost unbiased (see W2..1\i . each node is a 
random snapshot of the ULIRG population. In this sec- 
tion, we explore the conclusions that can be drawn from 
this diagram and its coded variants. In so doing, we as- 
sume that the IRS spectra are a product, on average, of 
the power sources governing the total IR emission. This 
assumption is reasonable for the population as a whole, 
but may break down for individual objects. 

5.1. The Network 

The first conclusion we draw is that, while graph- 
theory based approaches are a powerful tool for visual- 
izing complex and heterogeneous datasets, they require 
a large number of nodes. Our study has 102 nodes, and 
yet the conclusions we can draw from Figure [3] alone are 
limited to the subdivision of the sample into two or three 
groups. Solely from Figure [3] we cannot say if groups A, 
B and C are phases in time, or phases in some other 
variable. Even if we assume they are temporal phases. 
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then it is impossible to say what the time-ordering of the 
groups is. If we increased the number of nodes by an 
order of magnitude, then the resulting increase in reso- 
lution of the network may lead to more insight, but we 
cannot test this hypothesis here. 

Turning to the coded variants of Figure [3l it is clear 
that we still lack a complete and homogeneous dataset 
for z < 0.4 ULIRGs. Several nodes lack morphological 
and/or optical spectral classifications, and most do not 
have a dynamical estimate of central black hole mass. 
This omission is serious, given that local ULIRGs are 
the most easily accessible templates we have for under- 
standing the high redshift ULIRG population. 

We now examine the implications from Figure [3] and its 
coded variants for the ULIRG population. First, we ex- 
amine the drivers behind the division of the network into 
groups A, B and C. Three important lines of evidence are; 
(a) the presence of nearly all the widely and moderately 
separated systems in group A, (b) the fraction of single 
nucleus systems in groups B/C is much higher than in 



group A, and (c) the generally lower black hole masses 
in group A compared to group B, though the number 
of black hole mass measurements is too small for this to 
have much weight. We therefore propose that groups A 
and B/C represent temporal phases and are not signifi- 
cantly determined by other factors. 

Next, we examine the likelihood of group C being a 
separate entity to groups A and B. Group C is unlikely to 
be outliers to group A as its nuclear separations and PAH 
strengths are different from the outliers on the left hand 
side of group A. Instead, it is similar to the first 'half 
of group B. Based on the lack of connections between 
groups B and C, we tentatively propose that group C is a 
distinct stage from groups A and B, but stress that this 
is not robust. For example, it is conceivable that a source 
could move from group A to group C and then 'jump' to 
group B, although the lack of a bridge node connecting 
C to B means such a jump phase is likely short; as we 
have ^ 100 nodes, the lack of a bridge node suggests a 
jump phase length of order 2% or less of the total ULIRG 
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Fig. 13. — The degree distribution of Figure[3l i.e. the histogram 
of number of nodes with a given number of connections, k. The 
hght grey line is a (scaled) Poisson distribution with a mean of 7. 
Our degree distribution is clearly not Poissonian, and so does not 
follow the degree distribution of a random graph. 



lifetime^*. 

Overall, we propose that groups A, B. and C are dis- 
tinct but overlapping evolutionary phases, with A occur- 
ring first, followed by B and/or C. If a merger remains 
a ULIRG for most of the duration of the merger, then 
we can also estimate timescales based on the number of 
objects in each group; phase A lasts just over half the 
lifetime of a ULIRG, and phases B and C last approxi- 
mately half and one third the duration of group A, re- 
spectively. We do not claim that a ULIRG starts at the 
left end of A and goes gradually to the right. Instead, 
we propose that ULIRGs start in group A, with posi- 
tion and intragroup movement determined by unknown 
factors, and then proceed to B or C. We also note that 
the tests of reliability in involving randomization of 
the TZ values and varying the extinction still gave a clear 
subdivision of the network into groups A and B/C, but 
groups B and C were somewhat 'blended' and less dis- 
tinct. We conclude that the division of the network into 
'A vs. B/C' is more robust than the division 'A vs. B 
vs. C'. 

If the network structure is governed by temporal evolu- 
tion, we can use the purely network based metric of Be- 
tweenness Centrality to make testable predictions. The 
Betweenness Centrality (which we term B) of a node 
is the number of shortest paths betw een other pairs 
of nodes that pass through that node (jFreemanI Il979l : 
lBrandesl[200l . A low B means the node is inconse- 
quential, while a high B means the node is an impor- 
tant junction. A suitable analogy would be airports; a 
regional airport would have a low B, while an interna- 
tional hub would have a high B. It is therefore plau- 
sible that a node in the network with a high B is an 
archetype of a 'transitional' phase that many ULIRGs 
pass through. Most of the nodes in Figure [3] have B val- 
ues in the range 100 < B < 400. Fifteen nodes have B 
values in the range 400 < B < 800, while the four un- 

If we assume a total ULIRG lifetime of 10* years then the 
jump phase would be < 2 X 10® years long. This is short but 
feasible; for example, some Wolf-Rayet stars are expected to live 
approximately this long. 



connected nodes have B — 0. Six nodes however have 
what appear to be unusually high B values; Mrk 231 
{B = 1850), IRAS 00275-2859 (1620), IRAS 03538-6432 
(1420), IRAS 05189-2524 (1260), and IRAS 14348-1447 
and Mrk 273 (both --1000). We propose that these six 
objects are examples of key evolutionary phases. Based 
on their positions in Figure [31 we speculate that Mrk 273 
and IRAS 14348-1447 are templates of group A objects, 
IRAS 03538-6432 is a template of an object transitioning 
from phase A to phase B, IRAS 05189-2524 is a classic 
example of an object on the boundary between groups 
B and C, and Mrk 231 and IRAS 00275-2859 are prime 
early B to late B type objects. 

5.2. The Groups 

We now turn to the properties of groups A, B and C. 
The duration of group A is hard to quantify, but a rea- 
sonable estimate would be from the start of the merger 
to around the time the progenitor galaxies physically co- 
alesce. This is based on the broad range of projected nu- 
clear separations in this group, from widely separated to 
single nucleus. The almost universally prominent PAH 
features suggest the IR emission is powered mainly by 
star formation, though this does not preclude the pres- 
ence of a luminous AGN. 

Group A is also highly interconnected, and, as de- 
scribed in is assortative. In other fields where as- 
sortative networks are seen, neighboring nodes tend to 
identify as common members of a group, and /o r have 
similar properties (see lLusseau fc NewmanI ([200l for an 
interesting example). We therefore propose that star- 
bursts in group A are similar, at least to the extent to 
not give rise to significant differences in the IRS spec- 
tra. We propose that outliers to this group are instead 
caused by heavy intrinsic obscuration (see for example 
§5.4p . and speculate that there are no large variations in 
stellar IMF or metaUicity from ULIRG to ULIRG. 

Phase B follows and overlaps with phase A. Based on 
the fact that nearly all the objects on the right hand 
side of group B have single nuclei, group B likely ends 
some time after the nuclei of the progenitors have co- 
alesced. We see three interesting trends as we move 
from left to right in this group. First, the PAHs de- 
cline in prominence, becoming negligible (with respect 
to the continuum) by the right hand side. Second, sil- 
icate absorption varies from strong to weak on the left 
side, but is universally weak or in emission on the right 
side. Third, the optical spectral types are varied in the 
first half but almost universally Seyferts in the second 
half. We therefore propose that the relative contribution 
from star formation to the mid-IR emission declines as we 
move from left to right, while the contribution from an 
AGN increases^^, until some objects at the end of this 
phase briefly become optical QSOs. This is consistent 
with previous studies which show that the AGN fraction 
increases w ith increasing m erger age (e.g. IVeilleux et al.l 
I2002ll2006l though see also lRigopoulou et al.lll999f ). and 
suggests either that the accretion rate has increased upon 
moving into group B, and/or that the central black hole 
has reached a 'threshold' mass for luminous AGN activ- 
ity. We further propose that the heterogeneity of phase 

In contrast to phase A, where we make no claims relating a 
ULIRGs intragroup position to its evolutionary stage 
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Fig. 14. — Tests of robustness. The top panel shows the network diagram obtained if all 5151 TZ values are randomly changed by ±50. 
The bottom panel shows the network diagram obtained if the Ti values arc computed with a relaxed limit on the maximum foreground 
extinction of A9.7^m — 80, instead of Ag.r^m — 50. 




Fig. 15. — An example of a 'raw' Network diagram. This fig- 
ure was made using the same procedures as for Figure |3] but the 
positions of the nodes have not been subsequently adjusted. The 
nodes are coded by 11.2fim Equivalent Width, as in Figure [8l 



B arises from two factors. First is varying amounts 
of gas/dust driven into the nuclear regions. Second is 
AGN feedback; a luminous AGN can generate nuclear or 
galactic-scale winds, and the effects of these winds will 
vary substantially from case to case. 

Two examples lend support to our proposal that AGN 
become more luminous and less obscured as we move 
through phase B. First, IRAS 19254-7245 (object 81, also 
known as the Superantena) is located where we expect 
the AGN to be intrinsically luminous but still deeply 
buried, an expectation that app ears to be bor ne out by 
recent observations (Br aito et al.1120091 ). Second, 

Mrk 231 (object 8) is located where an AGN-driven wind 
may be expected, and indeed this object is thought to 
contain a starburst, an en ergetic AGN, and a nuclear 
outflow (jLipari et al.ll2005f ). 

Assuming that phase C is a separate evolutionary 
stage, then it is difficult to interpret as it contains a small 
number of objects. It seems to have similar properties to 
the first half of phase B, except perhaps for the nuclear 
separations, which may be smaller, on average, in phase 
C. Therefore, we propose that this phase also follows 
phase A, and that it is characterized mainly by waning 
star formation. We do not see evidence for a substantial 
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Fig. 16. — Schematic of the evolutionary scheme described in il5.ll 



increase in AGN activity in this group, and so propose 
that this phase is shorter in time than group B, and that 
its members are unhkely to become optical QSOs. 

We show a cartoon type diagram of this scheme in 
Figure [HI 

If all ULIRGs start off in group A, then the obvious 
question is what determines if they go to phase B or 
phase C?^" Broadly, there are two possible drivers; the 
dynamics/morphology of gas and dust (e.g. how much 
is available, and how efficiently it is channeled into the 
nuclear regions) , or seed black hole mass (models suggest 
the minimum seed bl ack hole mass for AG N activity in 
ULIRGs is - IO^Mq (iTaniguchi et aLinjOl l. If the lat- 
ter is the driver, and t he end-product of a ULIRG is an 
elliptic al galaxiy (e.g. iGenzel et al.ll200it iDasvra et al.l 
l2006bD then the antecedents of phase C should have 
smaller mass bulges than phase B. As there is no plau- 
sible evidence for a bimodality in BH/bulge masses in 
ellipticals, we think it more likely that the end products 
of B and C arc similar, except that some aspect of the 
merger dynamics of objects in B allows some of them to 
go through a brief optical QSO phase. 

This evolutionary picture fits well with recent studies 
of IR- luminous galaxies. The interconnected nature of 
Figure [3] implies that the starburst and AGN activity 
arises from a common physical mechanism, which tal- 
lies with imaging st udies of ULIRGs, which sh ow them 
all to be me rgers ( Surace et al.l l2000t ICui et al. 200ll: 
Farrah et aU 120011 : iBushouse etal l2002l : IVeilleux et all 
20G2„ .200611 The location of aU the QSOs in the 
sample at the end of group B, and their low B val- 
ues, suggests that few ULIRGs pass through a phase 
where they are simultaneously ULIRGs and Quasars, 
and/or that the ULIR G-QSO phase is brief, in agree- 
ment with recent work (iFarrah et al.ll200ltlTacconi et al.l 
I2002t iKawakatu et all l2006l . l2007t ) Our picture takes 
the idea that IR-luminous starbursts are present in most 
ULIRGs, w hile IR-luminous AGN are present in just 
under half (jGenzel et all Il998l : iRigopoulou et al.l Il999l : 

Assuming they cannot do both, but see il5.ll 



Tran et all 120011: iKlaas et all 120011: [Farrah et al] l2003l: 



Franceschini et al.l 120031 : IVega et all 12000) and extends 
it by providing (1) a single diagrammatic representation 
of the ULIRG evolutionary plane, (2) groupings into evo- 
lutionary phases, and (3) descriptions of the properties 
of these phases, including homogeneity, timescales, and 
power source. 

Finally, we note a peculiar aspect of the diagrams in 
Figure |5J the remarkable homogeneity of the 11.2/^m 
PAH strengths in group A, and the smoother gradient 
of 11.2/im PAH strengths through groups B and C, in 
comparison to the 6.2/^m PAH strengths. We do not 
have a plausible explanation for this. It could for ex- 
ample be highlighting an important part of the way in 
which PAHs diagnose star formation rates, or a subtle 
systematic error in our calculations. We do not consider 
this point further here, but highlight it as an interesting 
avenue to pursue in future work. 




Fork' Class 



Fig. 17. — Network d iagram, with the nodes color-coded by 
their 'Fork' classification l lSpoon et al.|[20071) . A copy of the Fork 
diagram has been embedded for reference. 



5.3. Comparison to other Mid-IR classification schemes 

As our evolutionary framework is based mostly on mid- 
IR spectroscopy, it is interesting to compare it to pre- 
vious work i n this field. A recent example is that of 
ISpoon et al.l ()2007f ) who published a mid-infrared based 
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Intrinsic spectrum A 



Intrinsic spectrum B 




A: fbi 
B: (l-f)bi 

(single relative scaling for all bins) 



A: fi 
B:(l-fi)bi 

(different relative scaling for each bin) 



Add cold dust 



Add Gaussian 
fluctuations 



P(A,B I same) 



Fig. 18. — Flowchart of the method used to calculate Ti. 
P ( A, B I different), and once to calculate P(A, B|same). 




For each pair of spectra, the flowchart is followed twice; once to calculate 



classification scheme (the 'Fork' diagram, their figure 1) 
for IR-luminous galaxies. The 'Fork' diagram and our 
network diagram overlap significantly in information use, 
as the Fork diagram uses both the silicate feature and the 
6.2^m PAH feature, but our diagram includes every other 
part of the spectrum, and is constructed using a funda- 
mentally different methodology. Our intention though is 
to compare the predictions from the two schemes, not to 
perform an independent check. 

In Figure [T7] wc code each node by its classification in 
the Fork diagram. We see a clear delineation. With a 
single exception, all the objects in group A reside in the 
upper branch of the Fork diagram (classes 2A/2B/2C, 
hereafter. Fork classes are given in italics). Group B on 
the other hand contains all the class 2A objects, nearly 
all the class lA objects, and a few 1B/2B/3A objects. 
Finally, group C contains the remaining class 2A''s, some 
iiJ's, and a few 2B''s and SA's. 

This indi cates that the two schemes are crudely in 
agreement. iSpoon et al.l (|2007[ ) suggest an evolutionary 
picture in which sources move up the diagonal branch in 
their figure I {2C - > 2B - > SB ~ > 3A), before either 



dropping vertically downwards {3 A — > lA) or via the 
slanted branch back to IC and on from there to lA/lB. 
In Figure [9] wc see a trend where the 2C/2B/3B sources 
lie on the left hand side, with the lA sources lying on 
the right hand side. Figure flTl however provides a more 
nuanccd diagnostic. From it, we can discern two distinct 
evolutionary paths, and that the IB and 2A classes are 
likely starburst/AGN transition classes, rather than just 
the 2A class. 

5.4. Notable Objects 

The positions of some famous objects are interesting. 
We describe some of them in this section. 

Arp 220 (object 7) is frequently used as a template for 
high redshift objects. Its presence in group A marks it as 
young. It is connected by only 4 edges, making it an atyp- 
ical example of a local ULIRG. This does not preclude its 
use as a template, but suggests it is not suitable to use if 
one is interested in determining if high redshift ULIRGs 
resemble local ones. From ! j5.11 its outlier status arises 
due to heavy intrinsic obscuration, which agrees with the 
properties of its mid-IR spectra; it shows strong PAHs, 
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Fig. 19. — Distributions of TZ calculated on spectral pairs with 
the same intrinsic spectrum but different dust extinctions. The 
black histogram is for those randomly generated spectral pairs 
with 5% Gaussian uncertainties while the grey histogram is for 
those pairs with 2% uncertainties. As expected, those pairs with 
"true" uncertainties less than those used in the calculation of 
are systematically shifted to smaller values. 
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Fig. 20. — The distribution of values of be resulting from maxi- 
mizing Eqn. IA30I for all possible pairs of spectra, fee is defined as 
the best estimate of the sum of the column densities of the two 
ULIRGs that obtains the best fitting intrinsic spectra. 



but also strong silicate absorption, and a steeply rising 
continuum at > 10/ini. 

Next, we consider IRAS F00183-7111 (object 11). This 
object is atypical (only two edges), and harbors an ob- 
scured AGN that is close to burning through its dust 
cocoon. An independent st udy comes to similar con- 
clusions; ISpoon et al.l ()2009f ). using the high-resolution 
IRS spectra to look for outflow signatures in fine struc- 
ture lines (i.e. information that is not contained in the 
low-resolution spectra) show that this object contains an 
obscured nuclear outflow driven by an AGN. 

Finally, the four objects in Table [1] that are not plot- 
ted in Figurc[3](16,41,48,79, one could also include object 
2 here as its single connection depends on the assump- 
tions made when computing TV) are interesting because 
their lack of connection is unusual. We defer a com- 
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Fig. 21. — The dust attentuation factor, e~*'absW^ ^ for an ar- 
bitrarily chosen value of x. 

plete discussion on these sources to a future paper, and 
here only briefly describe some possibilities. Part of the 
reason for this may be luminosity dependent; object 16 
(IRAS 00397-1312) is the brightest in the sample, albeit 
by a small margin, and object 41 (IRAS 07145-2914) is 
the faintest by over an order of magnitude. It may also 
be due to a combination of unusual spectral properties; 
object 16 (IRAS 00397-1312) has the deepest CO absorp- 
tion feature of any object in the sample, and object 79 
(IRAS 18030+0705) has extraordinarily strong PAH fea- 
tures. It is possible therefore that these objects represent 
either a very brief and/or very rare evolutionary stage, 
or that the merger dynamics are highly atypical in some 
way. 

6. CONCLUSIONS 

We have taken a large mid-infrared spectroscopic 
database of low redshift ULIRGs, and applied to it two 
novel analysis methods; (1) a Bayesian based estimator 
of similarities between pairs of spectra that takes into 
account every spectral resolution clement, marginalizing 
over luminosity, foreground obscuration and instrument 
noise, and (2) a visualization algorithm based on force- 
directed networks that efficiently presents these similari- 
ties across the sample simultaneously. We combine these 
results with archival data to propose, with some reserve, 
the following evolutionary description for ULIRGs in the 
low-redsliift Universe: 

1 - The IR emission in ULIRGs is consistent with being 
driven by a single underlying physical process. We see 
no evidence for multiple, separate evolutionary tracks. 
There is however evidence for at least two and possibly 
three evolutionary sub-phases. 

2 - The first phase (phase A), through which all 
ULIRGs go, lasts from the initial encounter until approx- 
imately coalescence. The IR emission arises mainly from 
star formation, with a contribution from an AGN in some 
cases. The highly interconnected, assortative nature of 
phase A suggests that there is little variation in star- 
burst parameters from object to object, with observed 
variations in the IRS spectra instead caused largely by 
differing foreground obscuration. 

3 - At around the time the progenitors start to coalesce, 
a ULIRG can branch off into one of two phases. We 
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suggest that the track a ULIRG takes depends primarily 
on the initial impact parameters and dynamics of the 
merger and the availability of gas/dust, and to a lesser 
extent on the masses of the central black holes in the 
merger progenitors. 

4 - The first of these two phases (phase B) lasts approx- 
imately half the length of phase A. The relative contribu- 
tion from star formation to the mid-IR emission declines 
as we move from left to right, while the contribution 
from an AGN increases, until some objects near the end 
of this phase briefly become optical QSOs. Phase B is 
less interconnected and more heterogeneous than phase 
A, implying that more than just foreground obscuration 
is driving the shape of the IRS spectra. We propose 
that this increased heterogeneity arises from two factors; 
(1) varying amounts of gas/dust driven into the nuclear 
regions by differing merger dynamics, and (2) feedback 
effects from AGN-driven winds. 

5 - The second phase (C) lasts about one third the 
length of phase A. It is similar to phase B in that the 
mid-IR spectra are heterogeneous, but the decline in lu- 
minosity of the starburst relative to the AGN is less pro- 
nounced. Few or no systems on this track pass through 
a QSO phase. 

6 - We use the graph-theory based metric of node be- 
tweenness centrality to identify six ULIRGs that may be 
archetypes of key points in this evolutionary cycle. We 
propose that Mrk 273 and IRAS 14348-1447 are examples 
of phase A objects, IRAS 03538-6432 is a prime example 
of an object transitioning from phase A to phase B, IRAS 
05189-2524 is an example of an object on the boundary 
between phase B and phase C, and Mrk 231 and IRAS 
00275-2859 are prime early B to late B type objects. 

7 - The 11.2/im PAH feature appears to be a re- 



markably good diagnostic of the evolutionary phase of 
a ULIRG, more so than the 6.2/im PAH feature or the 
9.7/im silicate absorption feature. Even though we are 
using the entire spectral range to construct the network, 
all the objects in group A have prominent, homogeneous 
11.2//m Equivalent Widths, which then smoothly decline 
as we move left to right through groups B and C, until 
we reach the right hand side of both groups where the 
PAHs are negligible. 
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APPENDIX 
MATHEMATICAL DETAILS 

Introduction 

The formalism for determin i ng th e degree of similarity between two data vectors using a Bayes factor was first 
developed in detail in I Jeffrevi (|1961f ). who quantifies the belief that two flux measurements using different detectors 
with different acceptances are measuring the same flux. This, in essence, reduces to whether or not the differences 
between the two measured fluxes are what is expected if the assumptions about the acceptances of the detectors are 
correct. 

We here extend this formalism in two ways; (1) we extend the calculation to involve not just one measurement per 
experiment, but many measurements from many wavelength bins {i.e. a spectrum), (2) we incorporate the ability to 
marginalize over external variables (e.g. luminosity, dust extinction) so that contributions from those variables to the 
TZ values can be accounted for. The procedure is shown as a flow chart in Figure [T51 

Method 

We start by defining the flux in the i*'^ wavelength bin of spectrum A as //foi, where hi ranges from to oo, and 
from to 1. If the fluxes are given in photon counts, then hi is the mean number of photons in the i*'' bin of 
both spectra, // is the probability that the photons are emitted from source A, and 1 — is the probability that the 
photons are emitted from source B. The number of photons in the i*^ bin of spectrum B is therefore (1 — f^hi. An 
analogy would be collecting balls into two receptor bins with different probabilities of accepting an individual ball; in 
this case, hi is the mean number of balls that will enter both bins, f[ is the relative acceptance of one bin, and (1 — //) 
is the relative acceptance of the other. 

If the two spectra are identical, then we will see equal contributions from the two spectra to each bin: 



f[b, = (1 - f^h, (Al) 

and therefore it follows that: 
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where / and (1 — /) are the 'true' fluxes emitted by sources A and B, respectively. In other words, to ask whether or 
not the sources have the same intrinsic spectrum is equivalent to asking whether or not = /. 

Now, if the two spectra differ, but only by a multiplicative scaling factor, and if Na and Nb are the normalizations 
for spectra A and B, respectively, then: 

Note that / would be the same for all the bins if the two sources were the same (that is, if they differed only by a 
normalization factor). We will not assume a specific value of / a priori; hence it must be marginalized {i.e. integrated 
out, thereby accounting for all possibilities for /). It is important to note that if the sources are intrinsically different, 
then is not constrained by /. 

It is now straightforward to expand Eqn.[T]in terms of /, and bf. 

^_ P(A,B [different) 



P(A,B|same) 
lo df dbP{A,B\f', b)P{f', b|different) 
df /J" dfoP(A,B|/, b)P{f, 6|same) 



(A4) 



where 6|same) and 6 1 different) are the prior densities, which encode any information that might be known 

about b, /', and / before the data are taken. For a Bayes factor to be well-defined the priors that enter the calculation 
must be proper probability densities; that is, they must integrate to one. However, as noted below, for this problem 
we are able to finesse this restriction. 

Let us consider each of the terms in Eqn. IA4I in turn. If the counts in the i"^ bins of the two spectra follow Poisson 
statistics, then: 

P(A,B|/', 5) = n ^-^'^ , (A5) 

i—l 

where D is the number of wavelength bins in both spectra and f^bi and (1 — fi)bi arc interpreted as the mean photon 
count in the i'^ bin of spectrum A and B, respectively; and Ni and Mj are the numbers of photons that are present in 
this bin of spectra A and B, respectively. Note that // is allowed to vary independently in each bin. If the spectra are 

different, as is assumed to be the case when we are calculating P(A,B|/', 6), then there is no constraint on the relative 
normalizations of the two spectra from bin to bin. 
When calculating the probability of obtaining the two spectra given that they are the same, we must consider: 

P(A,B|/, 6) ^ n m/ • 

i—l 

In contrast to Eqn. lASi the spectra here are only expected to have different overall normalizations and fluctuate 
according to Poisson statistics. 

We now consider the priors in Eqn. IA4l where we shall assume flat priors for /j', / and 6, defined initially on compact 
sets. We also assume that and b are independent a priori as arc / and b. That is: 

P{f, fo|samc) = P{f\b, same)P(6|same) (A7) 

= P(/|same)P(&|same). (A8) 

We take the prior for / to be: 

P(/|same) = (A9) 

Jmax Jniin 

22 



where — > and fr,iax 1 • Similarly: 



P(6|same) = — -jj (AlO) 



We take limits of the relevant variables following the advice in IJavnes fc Bretthorsti 1)20051 ) . One reason for doing so is to handle 
variables which in the limit are defined on [0, oo). A flat prior is improper in this limit, that is, it does not integrate to one. 
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where bmin — * and bmax — > oo. And so: 



6|samc) = — ^. (All) 



1 

{frnax fvairi)(praax ^rnm) 

Similarly: 



P(/',&|diffcrcnt) = P(/'|Mifferent)P(&|same) (A12) 

= P(/'|different)P(6|same), (A13) 

and likewise for fl: 

P(/'|different) = — ^ (AM) 

\J max J miri/ 

where /4im ~^ fl^^^ax ~^ 1- Finally, wc assume: 

P(6|different) = — (A15) 

Therefore: 

P(/', 6|different) = _ .1 ^^^^^ 

\Jmax JminJ K'^rnax fJmm ) 

As noted above, in principle the priors should be proper. However, for this problem the normalization factors in the 
densities P(/, &|same) and P(/', 6|different) cancel in Eqn. lA4l and we can finesse the issue of improper priors. But this 
is only because fmax = f'max ^ 1 and /„i„ = f'^^^ , 

In general, one wants to be sure that in higher dimensional problems such as the one considered here, TZ behaves 
properly. For instance, suppose that the upper limit on / was not 1 but p and, similarly, the upper limit on f[ was 

^. The priors on h would cancel, but we would be left with an overall constant p/^^ that depends on the number of 
dimensions, D. There are ways to en sure that TZ behaves as would be expected. For instance, one can use a method 
proposed bv lBerger fc Pericchil (|2001[ ). where the priors are such that: 

j J dfdbPif, b\sa.iRc) = J J df'dbP{f,b\diScTcnt) = 1 (A17) 

and therefore TZ is guaranteed to make sense. Another is simply to check that TZ exhibits the behavior that would be 
expected from a Bayes factor. For instance, one would expect that as the uncertainties decrease {i.e. Ni oo and 
Mi — ^ oo), a larger fraction of pairs with different intrinsic spectra would have TZ < 1 and a larger fraction of pairs 
with similar intrinsic spectra would have TZ > 1. We will come back to this point later when wc show that the TZ used 
to calculate the similarity of spectra indeed behaves properly as the flux uncertainties go to O^''. 
Collecting the terms, the full Bayes factor for the Poisson case becomes: 

P(A,B|same) ^ ' 

(/,'6i)"-e'^'-'''- ((l-jOi.i)'''e-<^-A')'-i 



Ni'. Mi'. 



which can be evaluated to be: 



TZ- 



n 



i=l Ni+Mi + 1 



iD 



(A19) 



B{N +1,M +1) 

where B{n, m) is the beta function B{n, m) = r(ri)r(TO) /T{n + m), N = X^iLi a'^'^ -^^ ^ EiLi ^'^i- 

If Ti, were not a Bayes factor, this would not affect the results presented in this work as Ti. is only used as a (v ery good) discriminant 
between spectral pairs whose intrinsic spectra are different and those whose intrinsic spectra are the same (see Fig. llOt . 
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Gaussian Errors 

In many cases (including ours) the data are in units of flux, not photon counts, and the flux distribution for the i 
bin will follow a Gaussian distribution. In this case, the likelihood in Egn. lASl is re-expressed as: 



and the likelihood in Egn. lAGI as 



^ 1 [ftj-FAif 1 la-flbj-Fsif 



i=i 



where, for the i*^ wavelength bin, fbi and f^bi are the mean fluxes for source A; (1 — f)bi and (1 — f-)bi are the mean 
fluxes expected from sources B; Fai and Fsi are the measured fluxes; and aAi and asi are the errors on the measured 
fluxes, all in the i'^ bin. Plugging Ean. lA20l and IA21I into Egn. lAlSI 



7^ = 



UZ^IodflCdb. 



2-KaA i <^B 



" 



H^-J'i)l>i-FB 



Iq df U^=l ir dh 2^aL^B. 

where now bi can be integrated scmi-analytically (sec Scction lX?^ 



e ""i. e 



(A22) 



Including Extinction 

We also wish to account for the possib le presence of a screen of cold dust in front of our sources. To do this, we use 
the carbonaceous-silicate grain model of IWeingartner &: Draind ()2001f ). their model 'A', with Ry = 5.5 and the grain 
abundances per H increased by a factor of 1.42, although any dust extinction law can in principle be used. The effect 
of dust extinction is accounted for by a simple exponential factor, e"'""^^'^, where fc(A) is the absorption cross-section 
per mass of dust (err? j g) as a function of wavelength, and x is the column density [gj crr?\ 

Figure [211 shows the distribution of e"'""^^'^ for an arbitrarily chosen value for x. If xa and xb are the column densities 
for sources A and B, respectively, they are parametrized by xa = /e&e and xb = (1 ~ fe)be, where the parameters fe 
and be are analogous to the parameters // and bi for the flux distribution considered above. That is, fee is the sum 
of the dust columns in front of sources A and B, and fe is the fraction of the dust column in front of source A. The 
Bayes factor, Eqn. IA41 is then modifled to include these new parameters: 



^ ^ lo dfe Jo°° rfbe /o df dbP(A,B|/e, 6e, f,b)P{^, be, /' , g|different) 
df, db, df d6P(A,B|/e, fee, /, b) P{f,, b,,f, 6|same) ' 



and the likelihood of parameters fe, be, f and b is now: 



P{A,B\fe,be,f',b) = l[-= e -= e (A24) 

while the likelihood of parameters fe, be, f and b is; 



P{A,B\fe,be,f,b)^Yl-= e -= e (A25) 

f-J^ V2TTaAi v2TTaBi 

where ki is the absorption coefficient in the i*'* wavelength bin. In other words, the mean flux obtained from source 
A with no dust is f^bi and fbi for the 'different' and 'same' hypotheses, respectively. If dust is included^^, this flux 
would be attenuated by an additio nal fa ctor, e~'^'^=''=. 

In practice, to calculate TZ, Eqn. IA23l necds to be constrained as it has too many degrees of freedom. Let us assume 
that some maximum, fee max, has been chosen so that fee < fee max- The priors are then modified: 

•^'^ We have parametrized the dust in such a way that one could in principle ask the question whether or not two sources have the same 
column density. Under the 'same' hypothesis, fe = 1/2, while under the 'different' hypothesis fe would be marginalized. However, we do 
not consider this question here. 
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P(/e, /, 6|samc) - 6e) ^ ^^26) 



^rmn ) 



and: 



where 8 is defined as: 



P(/e, 6„ /', foldifferent) ^ ^^^r^ _ (A27) 

"emaxUmax JminI V'max 'Jmm) 

e(6emax-M = (n '"^'^ J (A28) 

l^'Ji f^e max ^ f^e 

Pntting everything together, the full likelihood then becomes: 

^ ^ /o d/e /o /o°° d6F(A,B|/e, b^, 6)P(/e, ^^e, /, 6|different) 

d/e /o" dh, /o rf/ /o /o" dbP{k,B\f,, 6e, /, 6)F(/e, 6„ /, 6|same) 

r: 

Jo #e Jo d&e n»=i Jo Jo exp i L j exp ^ ^^^-^^ L j 



/o /o'''' '^^^ /o^ # n^Li /o°° cf^^ exp (- '° "''""'o-ff ^'"'' ) exp 



[e-ki(l-f„)be (l-f)bi-FBi]^ 



Parameter hi can be integrated semi-analytically as shown in Section [A. 71 Both the numerator a nd the denominat or 
of Eqn. IA2"9l are calculated using adaptive Monte Carlo integration methods, VEGAS and miser (jPress "eFaIlfT99l) . 

There is a degeneracy in Ean. lA29] Suppose we are comparing two spectra which both have the shape of the 
extinction curve shown in Fig. 1211 In this case, the best estimate of the parameters can come from an 'intrinsic' 
spectrum b) that looks like the extinction curve and is not propagated through any dust, or that is flat and is 
propagated through dust (assuming that dust follows the extinction curve in Fig. l2ip . In practice, this situation is 

exceedingly unlikely, as the intrinsic spectra (which enter through 6), the relative normalizations (that enter through 
/) and the extinction correction all do the best they can to make a good fit. There are however examples where effects 
of the source distribution taking on characteristics of the dust to make the best fit can be seen; see e.g. the comparison 
between IRAS 12071-0444 and IRAS 17179+5444, where the extinction is close to zero, but the spectra follow the 
features of the extinction curve in Fig. [51] 

Ideally we would like a model for the intrinsic dust distribution in each ULIRG, such that all the dust, hot and cold, 
can be placed where it should be (in the extinction parameters fe and he) and the source where it should be (in / and 

h). Such a model does not however exist, and so the best we can do is attempt to account for differences solely due to 
cold foreground dust. While this is undoubtedly an improvement on not attempting to account for dust at all, it does 
mean that (without constraints on the shape of the source distribution) the contributions of dust and the source may 
not be proportioned correctly. 

Behavior of TZ 

As discussed previously, TZ must behave as a Bayes factor should - i.e. as OAi and (jBi 0, pairs with like 
intrinsic spectra should all have 7?.'s less than 1 and pairs with unlike intrinsic spectra should have TZ^s larger than 1. 
To this end. Fig. [TOl shows that the spectra with lower errors systematically have a larger separation in TZ than those 
with larger errors, indicating the likelihood is behaving as expected, although this separation is somewhat obscured 
by the decreasing accuracy /precision of the Monte-Carlo integration as the uncertainties decrease. 

To gather further intuition about the behavior of TZ, we study the case where the Gaussian (statistical) fluctuations 
in the data are smaller than what is assumed in the calculation. We generate pairs of spectra with the same intrinsic 
spectra but different dust extinctions. The two spectra are then fluctuated independently assuming 2% Gaussian 
uncertainties; TZ is then calculated for each pair assuming 5% uncertainties. As there would be less variation in the 
data than expected by the calculation, one would expect the spectra to look more similar than if the fluctuations in 
the data were, say, ~ 5%. Therefore, one would expect TZ to be less than if the fluctuations in the data were indeed 
~ 5%. FigdU shows that TZ behaves as we expect where those data with 2% uncertainties are systematically shifted 
to lower values than those with 5% uncertainties. 

Maximum A Posteriori (MAP) Probability 

We would like to maximize the a posteriori probability dcflned as the unmarginalized denominator in Eqn. IA29I 
(which accounts for the probability that the two spectra come from the sources with the same intrinsic spectrum). 
This is done for several reasons. First, maximizing this likelihood is one way of finding out how similar the two spectra 
are. Second, by maximizing a posteriori probability with respect to fe, be, f and b, one can obtain a best estimate 
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of the mean distributions of these parameters for spectra A and B. Third, by plotting the best estimate of b^, for all 
pairing of the spectra, we get an estimate for the upper limit for be- 

The third point is worth describing in more detail. The rationale for calculating the upper limit on be in this way 
is as follows. We start by (falsely) supposing that the intrinsic sources are the same, and that the only thing that 
differentiates them is an overall flux normalization and the amount of dust in front of them. We then calculate the best 
estimate for be by finding the MAP defined later in the text (Eqn. IA30P , thereby creating a distribution of the best 
estimates for be for our data sample. The maximum be found from this distribution is an estimate of the maximum be 
that might be found for any of the pairs of objects. If our assumption that the sources were the same was incorrect, 
the differences in the sources would randomize the be distribution to some degree, thus effectively only pushing the 
upper limit on be higher. That is, the "true" value of be is likely to be covered in the integration of be. 

One could pose two objections to this method; (1) for a given pair of objects, we are using it in the distribution of 
the best estimates for be to find the upper limit of be and then subsequently using that upper limit to calculate TZ 
with the very same pair; and (2) we are not using the shape of the distribution of be as a prior when calculating TZ. 
The first objection can be addressed by realizing that in our sample there are 5151 pairings, and the effect of using 
our knowledge of be from a given pair twice (once to develop the distribution of the best estimates for be, and once for 
calculating the result using our knowledge on the upper limit of be) is small. As for the latter objection, one can argue 
that the distribution is skewed by the fact that the intrinsic source spectra for the pairs are different and therefore a 
flat prior on be is at least as good an estimate of our a priori knowledge of be as a prior that takes into account the 
shape of the be distribution. 

Turning back to the a posteriori probability defined as: 



Maximizing Eqn. lA30l (or minimizing Eqn. IA3ip yields a distribution of be from all the pairs of objects considered 
in this work, from where one can empirically determine the limit on be, 6emax- The distribution is shown in Fig. 1201 
We chose 6emax = 0.02 as the nominal value as be < bemax contains the bulk of spectral pairs, and fee max = 0.03 to 
test the dependence of the network diagrams on d ust e xtinction (see Fig. fM)) . 

The column densities found by minimizing Eqn. l A3l1 should not be regarded as good estimates of the 'true' column 
densities because there is degeneracy in the intrinsic shapes of the spectra and the effects of extinction. For instance, 
two spectra that are intrinsically the same will yield the same measured spectra (within statistical uncertainties) 
provided there is no dust extinction. However, two ULIRGs may have different intrinsic spectra and suffer different 
amounts of dust extinctions, but the combination might conspire to make the measured spectra look similar. Also, 
as discussed above, ULIRGs with similar column densities might, in fact, render a low value for be while the intrinsic 
spectra take on characteristics of the dust as there are no constraints imposed on them. 



L(A,B|/e, be, f, b) = P(A,B|/e, 6e, /, h)P{fe,be, f, 6|same). 



(A30) 



Maximizing Eqn. l A30l is equivalent to minimizing: 




(A31) 



Integrations With Respect to b 



Integrating the Gaussians in Eqn. IA29] with respect to bi ([Gradshtevn fc Rhvzik|[2000l ). we obtain: 



22 



~^ ^( 1 ') 



^2 + ^-2 



2 ^ 1 



i-l , + ^ n 



fe-kifobcp (i_f)„-ki(l-fc)bcp 
2 ^ 2 



.f2„-2kifebe ^ (i_f)2 ^-2^ ( 1 -fp jb^ 2<t^ . 



/ ,/„-kjfcbc 



Fa. , (i-f/)o-''.(i-fo)bo 

^ ^! 



t.'2„-kifcbo (l_t.')2 
2 i 'T 



f2p-2k;feb„ (i_f)2 



/ tc-l'ifoboF. ; (i_f)„-ki(l-fo)bcFj,; 
7 ^ 7 



2^i^ . 



(A32) 



where erf (u) is defined as: 



erf(u) 



'T^ Jo 



(A33) 
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54 
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57 
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59 
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15 
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61 
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13 


24 
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+05 37 04.7 
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SI 


12.73 


43 


62 
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13 


36 
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SI 


12.47 


1,8,24,46,81,84 


63 


IRAS 
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13 


36 


50.7 
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7 


12.50 


14,32,37,40,46,85 


64 


IRAS 
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13 


47 


33.3 


+12 17 24.2 
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S2 


12.37 


1,18,57,75,97 


65 
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14 


09 


31.3 
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S2 


12.88 


9,60,71,73 


66 


IRAS 


14378-3651 


14 


40 


59.0 


-37 04 32.0 
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L 


12.07 


3,4,9,58,60,72,82,83,87,89,93,94 


67 


IRAS 


15001+1433 


15 


02 


31.9 


+14 21 35.1 
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S2 


12.48 


1,5,9,13,26,29,32,37,38,42,46,50,53,60,68,69,70,72,74,80,88,91,97 


68 


IRAS 
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15 


22 


38.0 


+33 31 35.9 
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H 
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69 


IRAS 


15462-0450 


15 


48 
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SI 


12.24 


1,8,32,37,46,67,68,88,97 


70 


IRAS 


16090-0139 


16 


11 
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0.134 


L 


12.58 
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71 


IRAS 


16300+1558 


16 


32 
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0.242 


L 


12.69 


9,20,38,47,51,65,70,73 


72 


IRAS 


16334+4630 


16 


34 
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L 


12.41 
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+35 49 01.7 


0.371 


L 


12.39 


25,36,46,47,65,70,71,85 


74 


IRAS 


17068+4027 


17 


08 


32.1 


+40 23 28.2 


0.179 


H 


12.33 


5,9,13,17,22,29,35,38,50,51,52,54,56,60,67,68.70,72.77,80.88 


75 


IRAS 


17179+5444 


17 


18 


54.2 


+54 41 47.3 


0.147 


S2 


12.30 


1,18,28,57,64,68,91,97 


76 


IRAS 


17208-0014 


17 


23 


22.0 


-00 17 00.9 


0.043 


H 


11.94 


3,33,44,47,83,87,89 


77 


IRAS 


F17252+3659 


17 


26 


57.8 


+36 56 39.5 


0.365 


7 


12.47 


22,50,74,91 


78 


IRAS 


17463+5806 


17 


47 


05.6 


+58 05 18.0 


0.309 


S2 


12.61 


3,4,42,53,60,70,72 


79 


IRAS 


18030+0705 


18 


05 


27.1 


+07 05 57.5 


0.146 


7 


12.16 


- 


80 


IRAS 


18443+7433 


18 


42 


54.8 


+74 36 21.0 


0.135 


7 


12.27 


12,17,29,32,35,38,51,52,54,56,67,70,74,88 


81 


IRAS 


19254-7245S 


19 


31 


21.6 


-72 39 22.0 


0.063 


S2 


12.19 


1,8,62,97 


82 


IRAS 


19297-0406 


19 


32 


21.3 


-03 59 56.3 


0.086 


H 


12.37 


3,4,9,38,58,60,66,72,83,86,87,89,92,93,94 


83 


IRAS 


19458+0944 


19 


48 


15.7 


+09 52 05.0 


0.100 


7 


12.34 


3,4,9,29,44,58,60,66,72,76,82,87,89,92 


84 


IRAS 


20037-1547 


20 


06 


31.7 


-15 39 08.0 


0.192 


SI 


12.52 


8,15,37,45,62,95 


85 


IRAS 


20087-0308 


20 


11 


23.9 


-02 59 50.7 


0.106 


L 


12.34 


9,25,36,42,44,51,63,72,73 


86 


IRAS 


20100-4156 


20 


13 


29.5 


-41 47 34.9 


0.130 


H 


12.52 


4,9,19,35,38,52,82,89,92 


87 


IRAS 


20414-1651 


20 


44 


18.2 


-16 40 16.2 


0.087 


H 


12.18 


3,4,31,33,66,76,82,83,89 


88 


IRAS 


20551-4250 


20 


58 


26.8 


-42 39 00.3 


0.043 


H 


12.00 


17,22,26,32,35,52,54,56,67,68,69,74,80,91 


89 


IRAS 


21272+2514 


21 


29 


29.4 


+25 27 50.0 


0.151 


S2 


12.10 


3,4,9,19,35,38,60,66,72,76,82,83,86,87,92,93 


90 


IRAS 


23060+0505 


23 


08 


33.9 


+05 21 29.9 


0.173 


S2 


12.55 


15,43,59,98,99,102 


91 


IRAS 


23128-5919 


23 


15 


46.8 


-59 03 15.6 


0.045 


H 


11.97 


1,17,18,22,50,67,68,75,77,88 


92 


IRAS 


23129+2548 


23 


15 


21.4 


+26 04 32.2 


0.179 


L 


12.43 


3,4,19,54,56,82,83,86,89 


93 


IRAS 


23230-6926 


23 


26 


03.6 


-69 10 18.8 


0.106 


L 


12.25 


3,4,6,21,66,82,89 


94 


IRAS 


23253-5415 


23 


28 


06.1 


-53 58 31.0 


0.130 


H 


12.37 


4,9,49,60,66,82 


95 


IRAS 


23498+2423 


23 


52 


26.0 


+24 40 16.7 


0.212 


S2 


12.51 


8,10,15,37,45,46,84 


96 


3C 273 


12 


29 


06.7 


+02 03 08.6 


0.158 


SI 


12.83 


100,102 


97 


Mrk 1014 


01 


59 


50.2 


+00 23 40.6 


0.163 


SI 


12.63 


1,64,67,68,69,75,81 


98 


Mrk 463 


13 


56 


02.9 


+18 22 19.1 


0.050 


S2 


11.80 


90,99,100,101,102 


99 


PG 1119+120 


11 


21 


47.1 


+11 44 18.3 


0.050 


SI 


11.29 


15,24,90,98,101,102 


100 


PG 1211+143 


12 


14 


17.7 


+14 03 12.6 


0.081 


SI 


11.76 


96,98,102 


101 


PG 1351+640 


13 


53 


15.8 


+63 45 45.4 


0.088 


SI 


11.88 


98,99 


102 


PG 2130+099 


21 


32 


27.8 


+10 08 19.5 


0.063 


SI 


11.60 


90,96,98,99,100 



^ Optical classification, taken mainly from IVeilleux et al.' figgg!), and also from'Bcrgvall ct al. 1986: Ar mus et al.l[T989t lAllcn ct al."1991': 'Mirabel ct al."1991': 'Sckiguchi fc WolstencrofH 
[19931 : ILeech et al.ll994l : IVeilleux et al.ll995l : IDuc et al.ll997.: .Lawrence et al.ll999i : .Stanford ct al. 2000, : Vcron-Cettv fc Veron.2001 : Kcwlcy ct al., 2001: Mcusingcr ct al. 2001. : Rupke et al.l 
l2005l : IDarline fc Giovanelli|[200^ ^ IR luminosities are either taken from lFarrah et all I I2003I) , or calculated from the IRAS fluxes using the same methods. Units are the logarithm of 
the rest-frame l-lOOO/xm luminosity, in Solar luminosities.'^ ID numbers of the ULIRGs that 'pair' with this object, i.e. that have log(7^) < O'^ Measured from the IRS spectrum. The 
redshift given bv [Stanford et aP ((2003) does not agree with the IRS data. 
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